home *** CD-ROM | disk | FTP | other *** search
/ Monster Media 1996 #15 / Monster Media Number 15 (Monster Media)(July 1996).ISO / prog_c / cuj0696.zip / DWYER.ZIP / QFLOAT / QSQRT.C < prev    next >
C/C++ Source or Header  |  1996-03-31  |  2KB  |  88 lines
  1. /* qsqrt.c        */
  2. /* square root check routine */
  3. /* full precision for 9 word mantissa */
  4. #include "qhead.h"
  5. #include <stdio.h>
  6. extern QELT qsqrt2[];
  7. /* constants for first approximation polynomial */
  8. #if WORDSIZE == 32
  9.  
  10. static QELT qsq2[NQ] =
  11. {-1,EXPONE-3,0,0xd14fc42f,0xe79ba800,0,0,0};
  12.  
  13. static QELT qsq1[NQ] =
  14. {0,EXPONE-1,0,0xe3e3c2ae,0x4c146700,0,0,0};
  15.  
  16. static QELT qsq0[NQ] =
  17. {0,EXPONE-2,0,0xa08bdc7d,0xd5ffe300,0,0,0};
  18.  
  19. #else
  20.  
  21. static short qsq2[NQ] =
  22. {-1,EXPONE-3,0,0150517,0142057,0163633,0124000,0,0,0,0,0,};
  23.  
  24. static short qsq1[NQ] =
  25. {0,EXPONE-1,0,0161743,0141256,046024,063400,0,0,0,0,0,};
  26.  
  27. static short qsq0[NQ] =
  28. {0,EXPONE-2,0,0120213,0156175,0152777,0161400,0,0,0,0,0,};
  29.  
  30. #endif
  31.  
  32. extern QELT qzero[];
  33.  
  34.  
  35. int qsqrt( x, y )
  36. QELT *x, *y;
  37. {
  38. QELT a[NQ], b[NQ];
  39. int i, m;
  40.  
  41. qmov( x, a );
  42. if( x[0] != 0 )
  43.     {
  44.     qclear(y);
  45.     if( qcmp( qzero, a ) == 0 )
  46.         y[0] = x[0];
  47.     else
  48.         printf( "qsqrt domain error\n" );
  49.     return 0;
  50.     }
  51. if( x[1] == 0 )
  52.     {
  53.     qclear(y);
  54.     return 0;
  55.     }
  56.  
  57. m = a[1];    /* save the exponent        */
  58. a[1] = EXPONE - 1;    /* change range to 0.5, 1 */
  59.  
  60. /* b = ( a * qsq2 + qsq1) * a + qsq0        */
  61. qmul( a, qsq2, b );    /* b = a * qsq2        */
  62. qadd( qsq1, b, b );    /* b += qsq1        */
  63. qmul( a, b, b );    /* b *= a;        */
  64. qadd( qsq0, b, b );    /* b += qsq0;        */
  65.  
  66. /* divide exponent by 2 */
  67. m -= EXPONE - 1;
  68. b[1] = (m / 2) + (EXPONE - 1);
  69.  
  70.  
  71. /* multiply by sqrt(2) if exponent odd */
  72. if( (m & 1) != 0 )
  73.     qmul( b, qsqrt2, b );
  74.  
  75.  
  76. /* Newton iterations */
  77.  
  78. for( i=0; i<8; i++ )
  79.     {
  80.     qdiv( b, x, a );
  81.     qadd( a, b, b );
  82.     b[1] -= 1;
  83.     }
  84.  
  85. qmov( b, y );
  86. return 0;
  87. }
  88.